Risk Management

Importing the Data

# import data
# returns are in %
# sum return percentages by day
# assumed equal investment weight in each firm
read_func <- function(value, I)
{
  filename = paste(value, ".csv", sep="")
  ret = read.csv(filename, header=TRUE)
  
  firms = ret[["PERMNO"]]
  count = length(unique(firms))
  port = count*I
   
  ret_sum = aggregate(.~DATE, data=ret, FUN=sum)
  ret_sum = ret_sum[order(as.Date(ret_sum$DATE, "%m/%d/%Y"), decreasing=FALSE),] 
  ret_date = unique(ret_sum[["DATE"]])
  
  keeps = c("RET","DATE")
  ret_sum = ret_sum[keeps]
  
  # gets average
  i = 1
  while(i < nrow(ret_sum))
  {
    ret_sum[i,"RET"] = ret_sum[i,"RET"]/count
    i = i + 1
  }
  
  data = list("portfolio" = port, "returns" = ret_sum)
  return(data)
}

invest = 1000000 # investment per firm
file1 = "returns_main"
file2 = "returns_comp"
file3 = "returns_main_hist"
file4 = "returns_comp_hist"

# reads files, gets list of returns and value
data_m = read_func(file1,invest)
data_c = read_func(file2,invest)
data_mh = read_func(file3,invest)
data_ch = read_func(file4,invest)

Histogram of Returns

# histogram function with normal dist overlay
vec_histogram <- function(x, names)
{
  h = hist(x, breaks=(length(x)/50), col="red", xlab=names$xlabel, 
          main=names$title) 
  xfit = seq(min(x),max(x),length=40) 
  yfit = dnorm(xfit, mean=mean(x), sd=sd(x)) 
  yfit = yfit*diff(h$mids[1:2])*length(x) 
  lines(xfit, yfit, col="blue", lwd=2)
}

# lists with labels for histograms
hist_label_m = list("title" = "Main Period Returns w/ Normal Curve",
                       "xlabel" = "Daily Returns")
hist_label_c = list("title" = "Comparison Period Returns w/ Normal Curve",
                       "xlabel" = "Daily Returns") 

# write histograms for historical returns
vec_histogram(data_m$returns[["RET"]], hist_label_m)

vec_histogram(data_c$returns[["RET"]], hist_label_c)

VaR, $VaR and Expected Shortfall

#function to calculate var
var_calc <- function(returns,port,a)
{
  r_vec = returns[["RET"]]
  vol = sd(r_vec)
  cumdist = qnorm(1-a,0,1)
  
  var = abs(quantile(r_vec,1-a))
  dol_var = abs(quantile(r_vec,1-a)*port)
  exp_short = vol*dnorm(cumdist)/(1-a)

  data = list("VaR"=var, "$VaR"=dol_var, 
              "Expected_Shortfall" = exp_short)
  return(data)
}

# confidence interval for var
conf = 0.95

# var calculations
var_m = var_calc(data_m$returns,data_m$portfolio,conf)
var_c = var_calc(data_c$returns,data_c$portfolio,conf)

# prints out var Variables
var_m
## $VaR
##        5% 
## 0.0102009 
## 
## $`$VaR`
##      5% 
## 2040179 
## 
## $Expected_Shortfall
## [1] 0.01635318
var_c
## $VaR
##          5% 
## 0.008870062 
## 
## $`$VaR`
##     5% 
## 283842 
## 
## $Expected_Shortfall
## [1] 0.01550245

Daily Models

# risk metrics one day modeling
oneday_f <- function(returns,hist_returns,a)
{
  # initial variance using historical data
  # variance based off of previous 10 days
  hist = hist_returns[["RET"]]
  variance_0 = var(hist[length(hist)-10:length(hist)])
  
  # variables for while loop
  lamda = 0.94
  i = 1
  variance = c(0)
  var = c(0)
  exp_short= c(0)
  cumdist = qnorm(1-a,0,1)
  
  while(i <= nrow(returns))
  {
    variance_1 = lamda*variance_0 + (1-lamda)*((returns[i,"RET"])^2)
    variance[i] = variance_1
    var[i] = -1*sqrt(variance_1)*cumdist
    exp_short[i] = sqrt(variance_1)*dnorm(cumdist)/(1-a)
    
    variance_0 = variance_1
    i = i + 1
  }
  returns$variance = variance
  returns$VaR = var
  returns$ExpShort = exp_short
  return(returns)
}

oneday_m = oneday_f(data_m$returns,data_mh$returns,conf)
oneday_c = oneday_f(data_c$returns,data_ch$returns,conf)
# garch model
garch_f <- function(returns,hist_returns,a)
{
  # initial variance using historical data
  # variance based off of previous 10 days
  hist = hist_returns[["RET"]]
  variance_0 = var(hist[length(hist)-10:length(hist)])
  
  # function to solve for garch model variables
  x.g = garchFit(~garch(1,1),returns[["RET"]])
  # summary(x.g)
  coef(x.g)
  
  # variables for loop and garch model
  i = 1
  variance = c(0)
  var = c(0)
  exp_short= c(0)
  cumdist = qnorm(1-a,0,1)
  alpha = coef(x.g)[3]
  beta = coef(x.g)[4]
  omega = coef(x.g)[2]
  
  while(i <= nrow(returns))
  {
    variance_1 = omega + beta*variance_0 + alpha*((returns[i,"RET"])^2)
    variance[i] = variance_1
    var[i] = -1*sqrt(variance_1)*cumdist
    exp_short[i] = sqrt(variance_1)*dnorm(cumdist)/(1-a)
    
    variance_0 = variance_1
    i = i + 1
  }
  returns$VaR = var
  returns$ExpShort = exp_short
  returns$variance = variance
  return(returns)
}
garch_m = garch_f(data_m$returns,data_mh$returns,conf)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          2767
##  Recursion Init:            mci
##  Series Scale:              0.007927996
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V    params includes
##     mu     -0.32110904   0.321109 0.0321109     TRUE
##     omega   0.00000100 100.000000 0.1000000     TRUE
##     alpha1  0.00000001   1.000000 0.1000000     TRUE
##     gamma1 -0.99999999   1.000000 0.1000000    FALSE
##     beta1   0.00000001   1.000000 0.8000000     TRUE
##     delta   0.00000000   2.000000 2.0000000    FALSE
##     skew    0.10000000  10.000000 1.0000000    FALSE
##     shape   1.00000000  10.000000 4.0000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1 
##      1      2      3      5 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     3748.5188: 0.0321109 0.100000 0.100000 0.800000
##   1:     3735.3607: 0.0321103 0.122741 0.0967439 0.811486
##   2:     3726.3122: 0.0321096 0.121029 0.0739978 0.799681
##   3:     3725.3035: 0.0321096 0.126442 0.0752300 0.805061
##   4:     3724.6123: 0.0321098 0.121272 0.0712490 0.809206
##   5:     3724.4035: 0.0321100 0.115851 0.0732870 0.814325
##   6:     3723.9937: 0.0321098 0.114808 0.0679513 0.819821
##   7:     3723.9352: 0.0321098 0.113600 0.0678421 0.819764
##   8:     3723.8892: 0.0321098 0.113488 0.0682050 0.820917
##   9:     3723.7844: 0.0321098 0.111549 0.0674358 0.822160
##  10:     3720.7091: 0.0321075 0.0789351 0.0513750 0.872125
##  11:     3719.6966: 0.0321056 0.0694367 0.0466286 0.886533
##  12:     3598.5937: 0.0318379 1.00000e-06 0.00453679 0.996889
##  13:     3571.3127: 0.0318379 1.00000e-06 0.00432796 0.996348
##  14:     3565.7068: 0.0318379 1.00000e-06 0.00488800 0.996199
##  15:     3563.4444: 0.0318025 1.00000e-06 0.00482748 0.995847
##  16:     3560.9257: 0.0317659 1.00000e-06 0.00501565 0.995948
##  17:     3559.5364: 0.0316923 1.00000e-06 0.00515500 0.995639
##  18:     3557.8900: 0.0315436 1.00000e-06 0.00533711 0.995742
##  19:     3557.1529: 0.0313948 1.00000e-06 0.00537041 0.995583
##  20:     3556.4992: 0.0310972 1.00000e-06 0.00552383 0.995613
##  21:     3555.9031: 0.0305018 1.00000e-06 0.00556246 0.995464
##  22:     3555.3510: 0.0293109 1.00000e-06 0.00572724 0.995467
##  23:     3554.8663: 0.0269292 1.00000e-06 0.00581758 0.995294
##  24:     3554.6053: 0.0221658 1.00000e-06 0.00606972 0.995246
##  25:     3554.5320: 0.0201772 1.00000e-06 0.00610699 0.995109
##  26:     3554.3335: 0.0221659 1.00000e-06 0.00610391 0.995166
##  27:     3554.1144: 0.0261431 1.00000e-06 0.00613606 0.995076
##  28:     3553.7914: 0.0340975 1.00000e-06 0.00609815 0.995168
##  29:     3553.5406: 0.0340975 1.00000e-06 0.00656361 0.994848
##  30:     3553.5246: 0.0340975 1.00000e-06 0.00656750 0.994868
##  31:     3553.5127: 0.0340975 1.00000e-06 0.00652820 0.994876
##  32:     3553.5037: 0.0341001 1.00000e-06 0.00653122 0.994890
##  33:     3553.4981: 0.0341051 1.00000e-06 0.00650697 0.994894
##  34:     3553.4936: 0.0341154 1.00000e-06 0.00650847 0.994904
##  35:     3553.4902: 0.0341360 1.00000e-06 0.00649107 0.994907
##  36:     3553.4871: 0.0341772 1.00000e-06 0.00649168 0.994915
##  37:     3553.4839: 0.0342596 1.00000e-06 0.00647828 0.994917
##  38:     3553.4800: 0.0344245 1.00000e-06 0.00647788 0.994923
##  39:     3553.4610: 0.0364572 1.00000e-06 0.00646208 0.994930
##  40:     3553.4572: 0.0370594 1.00000e-06 0.00643008 0.994953
##  41:     3553.4572: 0.0370690 1.00000e-06 0.00643197 0.994952
##  42:     3553.4572: 0.0370676 1.00000e-06 0.00643196 0.994952
## 
## Final Estimate of the Negative LLH:
##  LLH:  -9831.504    norm LLH:  -3.553128 
##           mu        omega       alpha1        beta1 
## 2.938719e-04 6.285312e-11 6.431960e-03 9.949521e-01 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                   mu         omega        alpha1         beta1
## mu     -8.092996e+07  1.343625e+11  9.167498e+05  1.977179e+06
## omega   1.343625e+11 -3.361234e+16 -3.473200e+12 -2.218284e+12
## alpha1  9.167498e+05 -3.473200e+12 -3.533850e+07 -4.289392e+07
## beta1   1.977179e+06 -2.218284e+12 -4.289392e+07 -7.180548e+07
## attr(,"time")
## Time difference of 0.03700209 secs
## 
## --- END OF TRACE ---
## Warning in sqrt(diag(fit$cvar)): NaNs produced
## 
## Time to Estimate Parameters:
##  Time difference of 0.2490139 secs
garch_c = garch_f(data_c$returns,data_ch$returns,conf)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          2781
##  Recursion Init:            mci
##  Series Scale:              0.007515565
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U           V     params includes
##     mu     -0.58499393   0.5849939 0.05849939     TRUE
##     omega   0.00000100 100.0000000 0.10000000     TRUE
##     alpha1  0.00000001   1.0000000 0.10000000     TRUE
##     gamma1 -0.99999999   1.0000000 0.10000000    FALSE
##     beta1   0.00000001   1.0000000 0.80000000     TRUE
##     delta   0.00000000   2.0000000 2.00000000    FALSE
##     skew    0.10000000  10.0000000 1.00000000    FALSE
##     shape   1.00000000  10.0000000 4.00000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1 
##      1      2      3      5 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     3675.5865: 0.0584994 0.100000 0.100000 0.800000
##   1:     3668.3849: 0.0584981 0.101288 0.113754 0.810783
##   2:     3659.5290: 0.0584967 0.0845767 0.118009 0.807664
##   3:     3600.9398: 0.0584844 0.0361826 0.174126 0.846706
##   4:     3582.7562: 0.0584806 0.0210711 0.188424 0.858273
##   5:     3545.9600: 0.0584691 1.00000e-06 0.178801 0.892129
##   6:     3543.6229: 0.0584666 1.00000e-06 0.173478 0.902708
##   7:     3534.9772: 0.0584640 1.00000e-06 0.161763 0.900971
##   8:     3482.0944: 0.0584437 1.00000e-06 0.0894399 0.938879
##   9:     3462.9763: 0.0584383 1.00000e-06 0.0684187 0.951615
##  10:     3443.9948: 0.0584276 1.00000e-06 0.0267274 0.977655
##  11:     3440.6693: 0.0584276 1.00000e-06 0.0275239 0.978485
##  12:     3440.1856: 0.0584274 1.00000e-06 0.0282446 0.977588
##  13:     3439.2256: 0.0584279 1.00000e-06 0.0303964 0.976776
##  14:     3438.7013: 0.0584283 1.00000e-06 0.0318794 0.975018
##  15:     3438.4746: 0.0584274 1.00000e-06 0.0340236 0.974185
##  16:     3438.4293: 0.0584199 1.00000e-06 0.0343723 0.973537
##  17:     3438.3831: 0.0583993 1.00000e-06 0.0338891 0.973985
##  18:     3438.3800: 0.0583566 1.00000e-06 0.0337647 0.974084
##  19:     3438.3796: 0.0583129 1.00000e-06 0.0337907 0.974093
##  20:     3438.3777: 0.0582692 1.00000e-06 0.0337435 0.974104
##  21:     3438.3717: 0.0578906 1.00000e-06 0.0334763 0.974299
##  22:     3438.3661: 0.0575114 1.00000e-06 0.0334895 0.974225
##  23:     3438.3550: 0.0571323 1.00000e-06 0.0335416 0.974252
##  24:     3438.3328: 0.0556157 1.00000e-06 0.0336053 0.974179
##  25:     3438.3223: 0.0540991 1.00000e-06 0.0336895 0.974151
##  26:     3438.3217: 0.0536471 1.00000e-06 0.0337230 0.974119
##  27:     3438.3217: 0.0536596 1.00000e-06 0.0337176 0.974125
##  28:     3438.3217: 0.0536630 1.00000e-06 0.0337186 0.974125
## 
## Final Estimate of the Negative LLH:
##  LLH:  -10162.93    norm LLH:  -3.654417 
##           mu        omega       alpha1        beta1 
## 4.033077e-04 5.648372e-11 3.371859e-02 9.741245e-01 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                   mu         omega        alpha1         beta1
## mu       -93462003.7 -2.458207e+09 -1.928012e+05 -2.336219e+05
## omega  -2458206558.4 -2.422838e+15 -4.340546e+10 -6.467782e+10
## alpha1     -192801.2 -4.340546e+10 -1.205340e+06 -1.660497e+06
## beta1      -233621.9 -6.467782e+10 -1.660497e+06 -2.499098e+06
## attr(,"time")
## Time difference of 0.03500199 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.368021 secs

Plotting of Daily Data

2000-2010 Data Period Variance

# graph variance, VaR for main data
time = c(1:nrow(garch_m))
year = substr(garch_m[1,"DATE"],7,10)
year = as.numeric(year)
time = time/252 + year
data = data.frame(garch_m,time)
data$rm_variance = oneday_m$variance
data$rm_VaR = oneday_m$VaR

plot_ly(data, x = ~time, y = ~variance, name = "Garch Model", 
        type = "scatter", mode = "lines", 
        line = list(color = 'rgb(205, 12, 24)', width = 1.5)) %>%
add_trace(y = ~rm_variance, name = "RiskMetrics Model", 
          line = list(color = 'rgb(22, 96, 167)', width = 1.5)) %>%
layout(title = "Main Data Daily Variance",
       xaxis = list(title = "Year"),
       yaxis = list(title = "Variance"))

2000-2010 Data Period VaR

plot_ly(data, x = ~time, y = ~VaR, name = "VaR from Garch", 
        type = "scatter", mode = "lines", 
        line = list(color = 'rgb(205, 12, 24)', width = 1.5)) %>%
add_trace(y = ~rm_VaR, name = "VaR from RiskMetrics", 
          line = list(color = 'rgb(22, 96, 167)', width = 1.5)) %>%
layout(title = "Main Data Daily VaR",
       xaxis = list(title = "Year"),
       yaxis = list(title = "VaR"))

1980-1990 Data Period Variance

# graph variance, VaR for comparison data
time = c(1:nrow(garch_c))
year = substr(garch_c[1,"DATE"],7,10)
year = as.numeric(year)
time = time/252 + year
data = data.frame(garch_c,time)
data$rm_variance = oneday_c$variance
data$rm_VaR = oneday_c$VaR

plot_ly(data, x = ~time, y = ~variance, name = "Garch Model", 
        type = "scatter", mode = "lines", 
        line = list(color = 'rgb(205, 12, 24)', width = 1.5)) %>%
  add_trace(y = ~rm_variance, name = "RiskMetrics Model", 
            line = list(color = 'rgb(22, 96, 167)', width = 1.5)) %>%
  layout(title = "Comparison Data Daily Variance",
         xaxis = list(title = "Year"),
         yaxis = list(title = "Variance"))

1980-1990 Data Period VaR

plot_ly(data, x = ~time, y = ~VaR, name = "VaR from Garch", 
        type = "scatter", mode = "lines", 
        line = list(color = 'rgb(205, 12, 24)', width = 1.5)) %>%
  add_trace(y = ~rm_VaR, name = "VaR from RiskMetrics", 
            line = list(color = 'rgb(22, 96, 167)', width = 1.5)) %>%
  layout(title = "Comparison Data Daily VaR",
         xaxis = list(title = "Year"),
         yaxis = list(title = "VaR"))